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Abstract. Recent galaxy redshift surveys have brought in a large amount 
of accurate cosmological data out to redshift 0.3, and future surveys are 
expected to achieve a high degree of completeness out to a redshift exceeding 
1. Consequently, a numerical programme for determining the metric of the 
universe from observational data will soon become practical; and thereby realise 
the ultimate application of Einstein's equations. Apart from detailing the cosmic 
geometry, this would allow us to verify and quantify homogeneity, rather than 
assuming it, as has been necessary up to now, and to do that on a metric level, 
and not merely at the mass distribution level. This paper is the beginning of a 
project aimed at such a numerical implementation. The primary observational 
data from our past light cone consists of galaxy redshifts, apparent luminosities, 
angular diameters and number densities, together with source evolution functions, 
absolute luminosities, true diameters and masses of sources. Here we start 
with the simplest case, that of spherical symmetry and a dust equation of 
state, and execute an algorithm that determines the unknown metric functions 
from this data. We discuss the challenges of turning the theoretical algorithm 
into a workable numerical procedure, particularly addressing the origin and the 
maximum in the area distance. Our numerical method is tested with several 
artificial data sets for homogeneous and inhomogcneous models, successfully 
reproducing the original models. This demonstrates the basic viability of such a 
scheme. Although current surveys don't have sufficient completeness or accuracy, 
we expect this situation to change in the near future, and in the meantime there 
are many refinements and generalisations to be added. 



PACS numbers: 98.80.-k, 98.65.-r 



1. Introduction 

In modern cosmology, many attempts have been made to determine the large- 
scale structure of the physical universe using constraints provided by cosmological 
observations and knowledge derived from local physical experiments. The most 
common approach is to adopt the postulate that the universe is spatially homogeneous 
on large scales — the Friedmann-Lemaitre-Robertson- Walker (FLRW) model. Hence 
using observational data to determine the few free parameters characteristic of such 
universe models has become the primary objective, and this overall framework has 
been presented in great detail in the literature. 

The cosmological principle (see [6] and [32]) expresses spatial homogeneity as a 
point of principle, whereas the copernican principle merely states that we are not 
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privileged observers. As pointed out in [9], the cosmological principle determines a 
complete universe model, but we cannot verify it fully due to the predictions it makes 
about parts of the universe far beyond our observations. Although the copernican 
principle only has implications for the observable universe, its validity can potentially 
be proven with observations. Despite this, it is the much stronger cosmological 
principle which is almost invariably assumed in practice. Certainly, a good argument 
for homogeneity is provided by the Ehlcr-Gcren-Sachs (EGS) theorem [7], which 
required exactly isotropic CMBR observations for all observers, and the 'almost 
EGS theorem', or Stoeger-Maartens-EUis (SME) theorem [30], which allows small 
anisotropics in the cosmic microwave background radiation (CMBR) and obtains an 
'almost FLRW geometry [21]. But our attempts to verify homogeneity should not 
stop there. 

The general assumption that the universe has a Robertson- Walker metric on 
the very large scale has served cosmology well, and is implicit in many calculations. 
But this makes verifying homogeneity rather tricky as there is a distinct danger of a 
circular argument. To analyse the cosmological data consistently requires the use of 
a non-homogeneous metric. 

Current and planned galaxy surveys are vastly increasing the amount of 
cosmological data available for analysis. Already in recent years there has been a 
dramatic improvement in the number of measured cosmological parameters and the 
accuracy of their values. Properties of the matter distribution have been well studied, 
but always with the assumption of a homogeneous background metric. As accurate 
cosmological data accumulates, the proper reduction and interpretation of the high 
redshift data will require knowledge of the cosmic geometry that is traversed by the 
light rays we observe. It will no longer be necessary to assume homogeneity, the data 
will make it possible to quantify the level of homogeneity on different scales. Hence, 
being able to prove the homogeneity of the observable region of the universe rather 
than assuming it in principle is a long term objective of the current project. There 
are of course a variety of methods for checking homogeneity, such as the Sunyaev- 
Zel'dovich effect, and it is important to pursue the full range of methods. 

It should be emphasised that radial homogeneity is far harder to prove clearly 
than isotropy. Our cosmological observations are restricted to our past null cone, 
which is a 3-dimensional slice through our 4-dimensional spacetime, and the expected 
variation of observations with redshift is affected by the cosmic equation of state, the 
evolution of the observed sources, and the geometry of spacetime. Disentangling these 
effects without assuming homogeneity is not a trivial exercise [13,22]. 

We wish to determine the spacetime geometry as far as possible from astronomical 
observations with minimal a priori assumptions. In principle a set of observations of 
the redshifts, angular diameters, and apparent luminosities of galaxies, as well as their 
number counts, combined with knowledge of the cosmic equation of state and the true 
diameters, luminosities and masses of the sources (and the evolution of these source 
properties), can be turned into metric information. The idea of reducing observed 
cosmological data to a metric was first explicitly discussed by Kristian and Sachs [17]; 
they examined how this could be done near our present spacetime position by deriving 
expressions in power series for some astronomical observations near the observer in 
a general metric, and demonstrated the difficulties faced in confirming homogeneity 
of the universe from observations. However, the problem of source evolution was 
barely addressed in their derivations. In the ideal observational cosmology program 
by Ellis and Stoeger and others [1-3,11,19,23,24,26-29], they took a shghtly different 
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approach to Kristian and Sachs as they aimed to determine what could and could 
not be decidable in cosmology on the basis of ideal astronomical observations, and 
so considered the limits of verification in cosmology. They worked with observational 
coordinates since all observational data are given, not on the usual spacelike surface of 
constant time, but rather on our past null cone, which is centred at our observational 
position on our worldline. Hence, the observational data can be used with ease in the 
implementation of any algorithm developed through using the Einstein field equations 
(EFEs) that are written in observational coordinates. 

Thus there has been a fair bit of theoretical work on how to determine the 
cosmic metric from standard obervations, but implementation has not been attempted 
and the two key issues of choosing appropriate numerical methods and handling real 
observational data have not been properly addressed. However Bishop and Haines [4] 
did make a numerical attack on the problem that was only partly successful. They 
treated the past null cone (PNC) as a time-reversed characteristic initial value problem 
(CIVP). Since the CIVP code is not intended to deal with a reconverging past null 
cone, their numerics blew up at the maximum in the diameter distance, and they 
were not able to extend past this poiniQ. As shown here, careful consideration of 
the nature of this maximum allows the integration to be continued to much higher 
redshifts. This is clearly a very big task, and will take years to develop into a rigorous 
algorithm generating believable results. We here describe the beginnings of such a 
procedure, necessarily simple at first. Our focus here is on turning the theoretical 
algorithm outlined in [22] into a workable numerical method, and thereby providing 
a demonstration of the viability of a key component of the problem. In tackling a 
problem of this magnitude, it is essential to start simply, which is the main reason 
for initially assuming spherical symmetry and zero cosmological constant. In the 
long term we envisage a much more general treatment. (We note however, that 
we are unavoidably at the centre of our past null cone, and spherical coordinates 
provide a natural description.) Even this simple first step provides some interesting 
challenges, which are discussed below. Our expectation is that the accuracy and 
especially the completeness of cosmological redshift surveys will be much enhanced in 
the coming years, and extracting the geometric implications of the observations will 
become possible. 



2. The Lemaitre-Tolman-Bondi Model and its null cone relations 

The general spherically symmetric metric for an irrotational dust matter source in 
synchronous comoving coordinates is the Lemaitre-Tolman-Bondi (LTB) [5, 18, 31] 
metric 

ds' = -de + f^Ym ^ ' 

where R'{t,r) = dR{t,r)/dr, and d^'^ = dO'^ +sm'^ 9d(p'^. The function R = R{t,r) 
is the areal radius, since the proper area of a sphere of coordinate radius r on a time 

I Their initial data on the PNC were the diameter distance, the 2-metric derived from image 
distortions, the matter 4-velocity derived from redshifts and proper motions, and the matter density 
derived from number counts. Knowledge of the true shapes, absolute luminosities and masses of the 
sources was assumed. They adapted an axially symmetric, zero pressure, zero rotation, zero A, CIVP 
code and tested it on spherically symmetric Einstoin-dc Sitter PNC data, finding it to be second 
order accurate. Although r and / are called luminosity distances throughout, it is evident from their 
equations and figure 4 that they are diameter distances. This work has not been followed up. 
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slice of constant t is AnR^. The function E ~ E{r) > —1/2 is an arbitrary function 
of the LTB model representing the local geometry. 

Solving the EFEs with A = gives us a generalised Friedmann equation for 
i?(i,r), 

and an expression for the density 

where M{r) is another arbitrary function of the LTB model that gives the gravitational 
mass within comoving radius r. Here E{r) also plays a dynamical role, it determines 
the local energy per unit mass of the dust particles. Equation ([2]) can be solved in 
terms of a parameter 77 = r](t,r), and a third arbitrary function ts(r) which is the 
time of the big bang locally: 



. 2 y V 6 y M 

R= (l-cosr;) , ,y-sin77=^ ^^^^^ ^; < , (6) 

for hyperbolic, parabolic and elliptic solutions respectivelj{|. Specification of the three 
arbitrary functions — M{r), E{r) and tBir) — fully determines the model. They 
constitute a radial coordinate choice, and two physical relationships. 



2.1. The past null cone (PNC) 

The notation and null cone solution used here were first developed in [20]. However, 
they chose to work with the parabolic LTB model, and hence, their gauge choice which 
locates the null cone of the observer at one instant of time is simpler. This gauge choice 
was generalised to all spatial sections, i.e. for all values of E, in [22]. In this latter 
paper, which we will call MHE, they gave a complete outline of the observer's null 
cone in the LTB model, and how one can relate the LTB model to observables using 
this more general gauge choice. Therefore, we follow the general outline given in MHE 
here. 

On the one hand specification of the three arbitrary functions is what determines 
the LTB model, and on the other the angular diameter distance and the redshift space 
number density are what is given on the observer's past null cone. Hence, we first 
need to locate the null cone, and then relate the LTB arbitrary functions to the given 
data. 

Human observations of the sky are essentially a single event on cosmological time 
scales, and as a result, being able to locate a single null cone is all we need here; no 
general solution is needed. On radial null geodesies, we have ds^ = = dO'^ = dcfP' . 

§ However, near the origin, it is the sign of RE/M rather than E that determines the type of sohition. 
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From ([1]), if the past null cone of the observation event {t — tQ,r 
is given hy t — t{r), then t{r) satisfies 



dt 



R'[i{r),r] 
VI + 2E 



dr — 



R' 



V1 + 2E 



dr 



0) (here and now) 



(7) 



We will denote a quantity evaluated on the observer's null cone, t = i{r), by a " ; for 
example R[t{r), r] = R, and we note that it is a function of r only instead of r and t. 
If we choose coordinate r in such a way that, on the past null cone of (to, r), we have 



i?' 



1 



(8) 



+ 

then the incoming radial null geodesies are given by 

i{r) ^to-r . (9) 

With our coordinate choice ([5]), the density ^ and the Friedmann equation ^ on 
the past null cone then become 

m' 



AirpR-" = 



VTT2E 



R 



2M 
R 



2E{r) . 



The gauge equation is then found from the total derivative of R on the null cone. 



dr dr 



and this, together with 
dR 



and pT|ij| . leads to 



dr 



-VTT2E= -R= - 



2M 



We can then obtain an expression for y/l + 2E{ 



W{r) = VTT2E 




2E{r) 




(10) 
(11) 

(12) 
(13) 



,(14) 



where a new variable W(r) is introduced. This expression tells us for which regions the 
spatial sections are hyperbolic 1 + 2E > 1, parabolic 1 + 2E = 1 or elliptic 1 + 2E < 1, 
based on data obtained from the null cone. We substitute (|14p into (fTO|) and rearrange 
it into the form 

2 



dM 

dr 



inpR 

dR 
dr 



M 



2npR^ 

dR 



dR 

dr 



1 



(15) 



The proper time from the bang surface to the past null cone along the particle 
worldlines is described by 



T(r) = i{r) -tsir) = to-r - ts 



(16) 



{I Although wc do not strictly know the sign of R, it is fairly safe to assume that it is positive on our 
past null cone on the large scales that we are considering. From now on we take the positive sign for 
the right hand side of equation mill . 
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2. 2. Redshift formula 

Since the cosmological observations are given in terms of redshift rather than the 
unobservable coordinate r, we need to express all the relevant quantities in terms of 
redshift z. In order to do this, the redshift formula is developed here. 
As shown in MHE and elsewhere, the redhsift in LTB models is 

ln(l + z)= f'"^ ^^^''^^ dr (17) 

for the central observer at r = 0, receiving signals from an emitter at r = rp„i- 

We need to find the redshift z explicitly in terms of r, R and p, which we will 
later relate to observables. We differentiate ([2]) with respect to r, then evaluate it on 
the observer's past null cone, and we find 



R' _ 1 



(18) 



RVTT2E i?2 

From we can get the derivative of W. Using equation PU)) to eliminate M' , and 
combining with equations (jl3[) and (|17p . it now follows that 

dz/dr I d'^R , , , ---^ 

- ' ' A^pR] —- (19) 




1 + z y dr2 

with z(r = 0) = 0. So theoretically we now have the redshift in terms of coordinate 
radius r from R(r) and p{r) directly if we integrate ()19|) with respect to r. 



2.3. The observables and the LTB model 

As explained, we are assuming a spherical metric with a central observer purely for 
purely pragmatic reasons — one does not tackle the full generality of a complicated 
problem all at once. For simplicity, we suppose there is only one type of cosmic 
source and wc only consider bolomctric luminosities as in MHE. Sec HcUaby [13] for 
a discussion of multiple source types and multicolour observations. It is assumed that 
the luminosity and the number density of each source can evolve with time; with 
the former written as an absolute bolomctric luminosity L, and the latter as a mass 
per source, ^. Isotropy about the Earth is assumed, and wc also assume that the 
universe is described by zero-pressure matter - 'dust', and galaxies or perhaps clusters 
of galaxies are taken as the particles of this dust. 

The two source evolution functions might naturally be expressed as functions of 
local proper time since the big bang, L{t) and h{t). However, one cannot be sure of 
the age of the objects at redshift z because the bang time is uncertain in a LTB model 
and also because the location of the null cone is uncertain. The proper time from bang 
to null cone will be a function of redshift, t{z), and the projections of the evolution 
functions on the null cone arc written as L and fi. Of course, r(z) is unknown until we 
have solved for the LTB model that fits the data. For the sake of simplicity, we will 
take L and fx to be given as function of z, and wc use / for the apparent luminosity 
and n for the number count observations. In practice, many observational studies of 
evolution express their results in terms of z. 

The area distance do (or equivalently the diameter distance d^) is the true linear 
extent of the source over the measured angular size, which is by definition the same 
as the areal radius of the source at the time of emission, i.e. R in the LTB model. 
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It multiplies the angular displacements to give proper distance tangentially and its 
projection onto the observer's null cone gives the quantity R. The luminosity distance 
is measurable if we know the true absolute luminosity of the source at the time of 
emission L. If the observed apparent luminosity is l{z), then from the reciprocity 
theorem [8] gives the relationship between the diameter distance cIa ~ R and the 
luninosity distance drl^. 



{l + zfR{z)^dL{z) 




(20) 



where c?io = 10 parsecs. 

Let the observed number density of sources in redshift space be n[z) per steradian 
per unit redshift interval, hence the number observed in a given redshift interval over 
the whole sky is 

Aim Sz . (21) 

Thus the total rest mass between z and z + Sz is 

Anfm 6z , (22) 

where fi{z) = I-i[t{z)] is the mean mass per source. Given the local proper density on 
the null cone p, the total rest mass between r and r + dr evaluated on the null cone is 



pd'^V — p — . dr 



(23) 



where d'^V is the proper volume on a constant time slice. Hence, from ([22|) . ([23|l and 
(0), we get 



dz 

R p ~ pn— . 

dr 



(24) 



2.4^. The differential equations 

Most of the equations developed above are given as differential equations (DEs), and 
numerically DEs are easy to work with. Therefore, we need a set of DEs. that will 
generate the values of r, M, E and hence <b from the observations. The LTB model 
implied by the observations is thus deduced. 

Transforming (|19p to be in terms of redshift z instead of coordinate r, we obtain 
the null Raychaudhury equation 

(l + z) + 



dR 



dz dr^ 



(l + z) + 



d^R 



dz^ 



dz 



) = -AnpRil + z) . (25) 



Substituting ^ into using the facts that "sr = 1 / "SF 

\3 



and 



d~r 



dr 
dz 



-jpr ■ rewriting it so that all terms involving R (and hence dR/dz and 
d^R/dz^) are on one side of the equation, we then have a second order DE for r{z): 

-I -1 



dz2 




(26) 



^ Note that equation (31) in MHE is incorrect. 
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Since we want to solve all our DEs (and hence get the values for our functions r, M 
and E) in parallel, we need to introduce a new variable such that we can rewrite (j26p 
as two first order DEs. 

We introduce a new variable 4> = 't'{z), defined by 



dr ^ 
dz ' 
and equation ((26)) then becomes 

# J f 



R 



dz \ 1 + z _dR_ 



(27) 



(28) 



dz 



If we transform (|14p into a function of z and take square root of both sides, using the 
inverse of ([77]) we then obtain 



dR (1 - ^ 



W = VI + 2E = + ^ . (29) 

2(j) 

dz 

We substitute ([M]) into ([TS]) and rewrite it as dM/dz instead of dM/dr, so together 
with equation (|^ we then get 

= 47r/tnVl + 2E . (30) 

dz 

Hence, equations p7 p -(|30 |) give us a set of coupled first order DEs that we use in 
order to generate the values for r{z), M{z) and W{z) (or E{z)) from the observational 
data. We can then obtain the values for fj, t and then the third arbitrary function ts (z) 
for the hyperbolic and elliptic cases by substituting these values into equations ([U, © 
and ([16]), with ([4]) and ([6]) evaluated on the null cone. However, there is a borderline 
case - the near-parabolic case, which is needed where the exact expressions become 
numerically intractable. See [Appendix A| for the details. Note that from equation 
P8|). if we know the values for R and Anjln, we can then solve for (j) independently 
without knowing the values of r, M and W] while solving for r, M and W depends 
on knowing (j). This property of the (j) equation will be very useful later on. 

2.5. Origin conditions 

At the origin of spherical coordinates, r = 0, we have R{t,0) ~ and R{t,0) ~ for 
all t. Hence, on the observer's past null cone equations (fT2|) and ([TS]) then become 

dR 
dr 

r=0 

and thus J? w r to lowest order near r = 0. From (fTUl) we then find that 

M' « Anpor^ , M « y ^Po^' ; (32) 
and from (|14p using a Taylor series for R, and working to second order in r, we get 



R' ^ VTT2E = 1 , (31) 



E 



ifd^RV 4 



(33) 



Obtaining the spacetime metric from cosmological observations 



9 



where cPR/dr^ is finite when r = 0. Note that equations (|32p and ([55)F1 give us 
A/cxS3/2. 

We can get the origin hmit for fin from equation (|24p . which is 

fin fa — ;^;-r — . (34) 



az g 

Also, if we substitute ([5T|) into (dH), after rearranging the expression, we then get the 
origin condition for dz/dr (and hence z) 

dz d'^R 



dr dz'^ 



(35) 



The behaviour of the DEs ([27| - ([30|) needs to be checked near the origin, i.e. 
2 — > 0. Since r and z have a hnear relation, we know that near the origin, R ^ z, 
AiTfln ^ z^, dR/dz is finite, d^R/dz^ = finite and AI ~ z^. Also we know that 
dz/dr{0) = l/{dR/dz{0)), so </> = finite. Hence, our DEs are well behaved as z ^ 0. 

2.6. Apparent horizon and the maximum in R 

In the early universe, the expansion is so rapid that the light rays that are headed 
towards us are actually getting further away. One can consider the set of photons 
that are all the same time away from observation to be an incoming wavefront. As 
the universe slows down, there comes a moment when the area of such a wavefront 
is stationary, and R{t, r) has reached its maximum value. The locus of such points 
for all incoming wavefronts is the apparent horizon. Hence, for the LTB model, the 
maximum of the areal radius (or diameter distance) down the null cone is where the 
null cone crosses the apparent horizon. We locate this point by the calculation below. 

Since the apparent horizon is the hypersurface in spacetime where R is 
momentarily constant, we put dR/dr = into and using ([5]), Q and pT|) . we get 



2M 
and hence 



+ 2E = VT+2E , (36) 



R = 2M . (37) 

We will see that this locus presents us with a particular difficulty in our numerical 
reduction of null cone data. Of course, in the case when the cosmological constant is 
not set to zero, and if we are considering both the future and the past horizon; the 
calculation and the analysis will be more complicated [12,16]. See also [14] for the 
observational significance of this locus. 

There are a few things worth considering here - our DEs become singular when 
we reach the maximum in the areal radius (diameter distance) R, i.e. dR/dz ~ 0. 
From equation (|37[) we know that at the maximum of R, we have R = 2M. If one 
looks at equations ([^5]) and (PO)) . it actually contains zero over zero at this point, and 
any numerical method will break down here. Further, from and (pS)) we can see 
that d^R/dz^ = —iirfln/R where dR/dz = 0. There is no problematic behaviour of 
4> here, as can be verified in the FLRW case. Hence, in order to carry our numerics 

+ Note that equation (22) given in MHE is incorrect. Their expression does not allow all three cases 
of E. 
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through the maximum of R, we need to perform a series expansion near this point for 
R{z), jiniz), <j>{z), M{z) and W{z), as given in (|B.ip - (|B.5p of [Appendix B| 

Here, we use Rmax to denote the maximum in R, and its corresponding z value 
is called Zm- The series are then written in powers of Az = z — Zm- From the R 
and (in data, we can easily determine the values of Rmax, fj-n{zm) and Zm, and thus 
the remaining i?(z) and fln{z) coefficients can be evaluated by simply performing a 
least squares fit using the data values near Zm- In order to obtain the expressions for 
the coefficients in the 0(z), M(z) and W{z) series, we need to substitute (jB.ip - (|B.5| 
from [Appendix B| into our DEs (HTj)-!!!!]). The detailed expressions can be found in 



[Appendix B[ 

From (jB.12[) - (|B.15[) we can see that all (/)(z) coefficients are determinable once 
we know the values of z„i and all coefficients of the R and Anjln fits. Using ([27|) and 
(|B.3p . the series expansion for r is simply 

r(z) = ro + cboAz + ^ Az^ + ^ Az^ + ^ Az^ + • • . , (38) 

where = r(zm) is the integration constant. 

Note that Mq is obtained directly from Rmax without any further information. 



The only problem is that expressions (jB.lOp . (|B.11|) and (|B.16p - (|B.18p in [Appendix B 
all depend linearly on Mi. Unfortunately, no information about Mi can be obtained 
when we carry out the series expansion, as one can see from (jB.Qp . Despite this, it is 
still possible to obtain a value for Mi by substituting a known value (from numerical 
integration), say Ma at z^, where z^ is some distance away from Zm, into (jB.4p as 
described in more detail below. 



3. Numerical procedure 

3.1. Data handling and numerical method 

The actual data we must use consists of redshift and apparent magnitude 
measurements for a large number of discrete sources, which must be sorted into 
redshift bins. In each bin we must calculate the total number of sources, 47rnJz, 
and the average value for the diameter distance, R. 

Now the above theory treats all physical and geometric quantities, such as the 
density and the metric, as continuous functions of position, while the available data is 
a discrete set of sources. Therefore it might at first seem one should fit a smooth curve 
to the data in order to proceed with the integration. However, numerical methods are 
not continuous cither, and any numerical method we might choose to solve our PNC 
equations with is based on discretisation of continuous DEs. Furthermore, we must be 
careful not to hide any inhomogencity by smoothing on too large a scale, or introduce 
unintended bias by inappropriate choice of smoothing function. We note that the 
process of calculating averages on the redshift bins already introduces a measure of 
smoothing and a basic smoothing scale. We also note that higher order integration 
methods that use data from several different z values will also have a smoothing effect, 
and this is an option we are keeping open. So although statistical fluctuations in the 
real data may need a further degree of smoothing, we prefer to keep it to a minimum, 
only introducing as much as necessary. We also argue that any extra smoothing that 
is required should be closely tied to the numerical integration scheme, and not merely 
ad hoc. 
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An important consideration in the choice of numerical method is that the right 
hand sides of the DEs ([77|) - ([5(H) contain not only the funtions being solved for, 
r{z), M{z) and E{z), but also the given data derived from observations, R{z) and 
A'KiJLn{z). Since the latter are only known at discrete z values (the mid points of the 
z bins), methods that allow adaptive step sizes are not appropriate, and similarly 
evaluations of 0(2:), r{z), M{z) and E{z) significantly above or below their correct 
values should be avoided because there is no way to find the corresponding values 
of the given data. A second consideration is that with real observations, there will 
be statistical fluctuations and measurement uncertainties, so there is a limit to how 
much improvement can be gained from using higher order methods. In line with our 
policy of not using a more complicated method than the situation demands, we found 
that, with bin size (= step size) Sz = 0.001, a second order Runge-Kutta method 
gave entirely satisfactory results when very accurate fake data was given for the data 
functions. These choices may change in the future once real data is used and as more 
factors are included. 

Once 4>{z), r{z), M{z) and E{z) are determined, then ?y(z), r(z) and tB^z) are 
easily obtained from the algebraic equations ([9|) and (|4|)-(l6]). However, at each discrete 
position, we are required to determine, numerically, which type of evolution we have: 
hyperbolic, elliptic or near-parabolic. Note that for the near-parabolic case, when E 
is small but not zero, we use equation (|A.6[) from [Appendix A\ As one might have 
noticed, equation (|A.6|) is in powers of R{2E)/M, and this factor can be evaluated at 
each discrete position since we know the values for R, E and M. Of course, the error 
for this approximation of the series expansion for r has to be small, say about lO^''; if 
we take (|A.6P up to order 3, this will give us ^^1^^ ~ 0.1. Hence, if ^^^f^ > 0.1, 

we use the hyperbohc case (|3]); if ^^^P < —0.1, we use the elliptic case and if 

-0.1 < -^^^ < 0.1 we use the near-parabolic case. 

A set of computer programmef0 were developed that generate the values for the 
LTB functions M, E and is. Below we give a brief summary of the order of the steps 
followed in our programmes. To obtain the mass, energy and bang time functions (M, 
E and Ib respectively) from observational data and source evolution, we proceed as 
follows: 

(i) take the discrete observed data for l{z, 9, 0) and n{z, 9, 0), divide it into redshift 
bins of chosen width 6z, and in each bin average or sum it over all angles and the 
bin width to obtain l(z) and n{z). We may wish first to correct the data for known 
distortions and selection effects due to proper motions, absorption, shot noise, image 
distortions, etc.; 

(ii) choose evolution functions L{z) and jl{z) based on whatever observations and 
theoretical arguments may be mustered; 

(iii) determine R{z) from L{z) and l{z) using (|20[) . this is then our first input 
data function and we have Anfln as our second input data function; 

(iv) numerically integrate the DEs ([27l) - (l30l) using the redshift bins as the basic 
step size, and the binned data for z, R and At: fin, thus obtaining r(z), M(z) and E{z); 

(v) solve for fj, T(r), and hence t_B(r) from (d])-® and ITBl) . with dH)-® evaluated 
on the null cone. Notice that L{t) and /x(r) could also be found in this step. 

However, at this early stage of development, steps (i)-(iii) use test data generated 
from a variety of model assumptions. In fact, step (iv) has four components, which 

* Matlab was used for our numerical work. 
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are summarised below: 

• Deduce the origin parameters and output for the first 3 data points, i.e. at z = 0, 
z = {l/2)5z and z = (3/2),5z. 

• Use numerical DE solvers - a second order Runge-Kutta method - for solving the 
DEs up to just before the maximum in R is reached. 

• Determine the point Za where the switch to the series expansion is made, evaluate 
all quantities to be matched such as Mq, Wq, etc, and extend the numerics through 
the maximum in R by calculating the series expansions of [Appendix B| for r, 0, 
M and E. 

• Evaluate another matching value for switching back from the series expansion to 
numerical integration, and continue to solve the DEs numerically up to the limit 
of the data, here set to z = 3. 

3.2. Data near the origin 

Although we have already discussed how the DEs behave near the origin, from any 
available cosmological data that we might use, there is no data available at the origin 
itself, and very little in the first few redshift bins. Therefore, a method of filling in 
this gap is necessary in order to provide the initial values needed by the numerical 
integration. 

If we average over all the data values within each bin, for example the R values 
within a given z bin, then the average R values that we use are located roughly in the 
middle of each bin. Thus we have the first value of at z = Jz/2. We can then get 
the discretised versions of dR/dz and d^R/dz^ from the first and second differences 
of R. It takes two 5z bins to get the first value of the first difference at z = dz and all 
the values are located at z = kSz for any positive integer k. However, it takes three Sz 
bins to get the first value of the second difference and hence the first value of d^R/dz'^ 
at z = 3(5z/2, therefore, all the values are located in the middle of each bin. 

Since we want to have a complete set of data at each z value, we take the average 
of the two neighbouring first difference data points to get all our data at the half 5z 
locations (in the middle of each bin). In doing so, we will not have data values at 
the origin or at z = 5z/2, since the first complete data set is at z = 3/2(5z. But we 
know that LTB is RW like near the origin due to the fact that it assumes spherical 
symmetry. For that reason, the series expansions of the RW expressions are used for 
finding the RW parameters that fit the data values at the origin and at z = ^z/2, 
i.e. we determine the central values of Ha (Hubble constant) and go (deceleration 
parameter) from the data near z = 0. So we take the standard RW expressions for 
R{z) and 47r /irt(z) given by equations (A.l) and (A. 2) in Appendix A of MHE, and we 
do series expansions of them near the origin, as detailed in [Appendix C[ The origin 
limits are r(0) = 0, (/)(0) = dR/dz{0) = 1/Ha, M(0) = 0, £^(0) = 0, R{0) = and 
d^R/dz^{0) = -i3 + qo)/Ho. 

3.3. Passing through the maximum in R 

As mentioned before, a series approximation is required in the vicinity of Rmax, where 
the DEs become singular. Here all 0(z) coefficients are determinable once we know 
the values of Zm and all coefficients of R and AiTjln. A set of (j) values can be generated 
from the series expansion by substituting a set of Az values into (jB.3p . for a z interval 
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that overlaps with our numerical results. Wc use the R and [in values of 180 redshift 
bins on either side of z.^ (this is a redshift interval of 0.361 which covers about 12% of 
the total redshift interval that we are considering here), and perform a least squares 
fit with these data to obtain all the coefficients for R and /in, and hence obtain all 
coefficients for 0. There is good agreement between the numerical and series values 
over a range of z values when plotted on the same graph, and there is one intersection 
point between the two curves before Zm (and also the closest to Zm)- The intersection 
points here are important since they are where we match values between the series 
expansion and the numerical integration for 0. This is what we choose for Zq, and the 
numerically derived M at Za becomes our Ma- 

Now that we know where Za is, we can get Mi from Za and Ma if we substitute 
them into equation (|B.4p . using Az = Za — Zm 

Ma = Mo + Mi{Za - z™) + M^iZa " Z^)^ + M^^Z^ " z^f + ■■■ .(39) 

where M2 and M3 are given by (|B.10[) and (jB.lip . Similarly, if we are matching the 
W values 

Wa=Wo + Wi [Za - Z„) + W2 [Za - Z^f ^ ■ ■ ■ , (40) 

where Wi and W2 are given by (jB.17p and (jB.lSp . Therefore, if wc match M then 



Ml 



Ma -Mo , 87r2(/i7l)o2 



R2 



Rmax 

(/in)ii?2 , 



{Za 



+ 



i?2 



{Za 



1 + Zm 

87r^(/<n)o(/t7i): 

Rrnax 
(A"-)2 



(A"-)o 



R3 



{Za 



4 
1 + 

i?2 



1 



[Za 

3(/i?i)o ' 3(/in)o(l + z,„) ' ° 

Alternatively, if W is used for the matching, we have 
i?2 , , , 47r(/in)o 



(/in)o 

n) 



(41) 



Ml = Wa 



(Za 



47r(/in)o 

3-R3 _ \2 



167r(/in)o 

3i?2 



(Za 



{fln)iR2 
167r(/in)Q^ 

,2 27r(/in)i 



(Za 



167r(/in)o(l + z,; 
1 



(Za ■ 



R2{Za - ZrnY 




A'K{fin)o 47r(/in)o(l + z„) 47r(/in)o^„ 



(42) 



All our functions should connect the numerical integration (z < Za) and series 
expansion (z > Za) parts at Za] and from Za, we can generate the corresponding 



ra and M^ easily. With all (f> coefficients known, a value for can be found from 
Using ((4T|) or (|42)) . a value for Mi can easily be determined. 

The purpose of doing a series expansion is to extend our numerics through Rmax , 
but once this is achieved, we need to switch back to numerical integration again. It is 
sensible if we connect at zj where z^ — Za = zj — z^^. 

Initially, we tried matching M at the two connecting points, Za and zj. However, 
after we compared r, (f>, M and W from our numerics with the correct curves generated 
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from the assumed model, the r, (p and M curves showed good agreement, but the W 
curve had jumps at the two connecting points. As anticipated, W is the least well- 
determined function. 

We then tried to match W at both connecting points. Although this removed 
the two jumps in the numerical W curve, it also reduced the accuracy of the W series 
expansion. A key consideration is that at Zm we actually know the value of M from 
(jB.SP if we know Rmax- In order to maximise the accuracy of our series expansion 
and minimise the jumps that appear in our W graph, we matched M at the first 
connecting point, and W at the second one. 

This approach does not leave any visible alteration in the M curve, the jump 
in W at the first connecting point is still present, but better accuracy for the series 
expansion is obtained and the second jump is avoided. One thing worth mentioning 
here is that we may need to shorten the z interval for the series expansion, since 
with inhomogeneous data, fluctuations will be present, so if the interval is too wide 
compared with the fluctuations, the accuracy for our series expansion will be lower. 
However, this problem will only be dealt with when it has shown a significant effect 
on the numerics. 

4. Testing the numerics 

In order to test our numerical procedure, we need fake "observational data" for which 
the correct results are known. Therefore we generated sets of "observational data" 
that would be produced in a selection of LTB universes. 

Although we did a full comparison of numerical output from our programme, 
M, E and ts, with the correct LTB functions for a variety of different models, both 
homogeneous and inhomogeneous, we cannot present all our results here. Therefore, 
we summarise the ones we did in Table [1] below and only present a complete set of 
plots from one homogeneous and one inhomogeneous model in the two subsections 
below. In order to avoid confusion, we call the Hq and qo used for generating fake 
data god and Hm', and the ones our numerical procedure extracts from the data qo 
and iJo from here on. Where the model is inhomogeneous, both pairs are the values 
at the origin, as explained in Section 3.2. 

^■l. Homogeneous models 

Amongst the homogeneous models we found that the near-parabolic models to have 
slightly lower accuracy. As noted above, we expect the output function with largest 
error to be = Vl 4- 2E. 

Below we present the results of the comparison for a homogeneous model with 
qod = 0.49 and iJod = 0.72. This is the case when we have a negatively curved 
universe, but very close to the flat case. In Figures [T][3l the curves plotted from our 
numerical output using the Runge-Kutta method and the ones from the generated 
RW data are in very good agreement with each other, and there is only 0.02106 % 
difference between the curves in Figure ID As can be seen in Figure |4l there is still a 
jump in W at Za although it is barely visible, while the jump at the second matching 
point is no longer visible. This justifies the earlier decisions to match first M and 
then W at the two connections between the numerical integration parts and the series 
expansion part. 
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Table 1. Summary of all the homogeneous and inhomogeneous cases used for 
full comparisons between our numerics and the generated fake data. 



Homogeneous cases 


Hyperbolic 


Near 


Elhptic 


(fc = +1) 


parabohc 


(fc = -l) 


qod ^ 0.45 


qad = 0.49 


qod = 0.8 


Hm = 0.72 


Had = 0.72 


ffod = 0.72 


qod = 0.1 


qad = 0.51 




Hm = 0.72 


Hnd = 0.72 




Inhomogeneous cases 


Hyperbohc 


Near 


Ehiptic 




parabohc 




Varying bang time 


Varying geometry or 


Strongly inhomogeneous 


with qod = 0.2 


energy with q^d = 0.52 


with qQd = 0.6 


Varying mass 






with gorf = 0.22 







The T curves in Figure [51 and the in Figure [SI show good agreement between 
the numerical output and the correct values. 



- 




LTB 
RW 



Figure 1. Results of r vs. z with Hq ^ 0.71999, 90 0.490004 and 5z = 0.001. 

The grey curve is the correct RW expression and the dotted black one is our 
numerical output using Runge-Kutta as the integration method. 
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LTB 
RW 



\ 



Figure 2. Results of vs. z with _ffo ^ 0.71999, go ~ 0.490004 and 5z = 0.001. 

The grey curve is the correct RW expression and the dotted blaclt one is our 
numerical output using Runge-Kutta as the integration method. 



LTB 
RW 



Figure 3. Results of M vs. z with Hq 0.71999, qo ^ 0.490004 and 5z = 0.001. 
Matching the M values at the first connection point and the W values at the 
second connection point. The solid grey curve is the correct RW expression and 
the dotted black one is our numerical output using Runge-Kutta as the integration 
method. 



Obtaining the spacetime metric from cosmological observations 



17 



RW 
LTB 



Figure 4. Results of W vs. z with k. 0.71999, (jo ~ 0.490004 and 5z = 0.001. 
Matcliing tlie M values at the first connection point and the W values at the 
second connection point. The solid curve is the correct RW expression and the 
dotted one is our numerical output using Runge-Kutta as the integration method. 



RW 
LTB 



Figure 5. Results of r vs. z with Ho « 0.71999, <jo ~ 0.490004 and 5z = 0.001. 
The solid grey curve is the correct RW expression and the dotted black one is our 
numerical output using Runge-Kutta as the integration method. 
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0.8 1- • RW 

LTB 

7 U Now 

0.6 - 
0.5 - 
0.4 - 



0.3 - 
0.2 - 
0.1 - 
. , 

-0.1 — 



Figure 6. Results of vs. z with witii Ho 0.71999, go ~ 0.490004 and 
Sz = 0.001. The thick solid grey curve is the correct RW expression and the 
dotted black one is our numerical output using Runge-Kutta as the integration 
method. The solid black line on the top is the current age of the universe. 



4-2. Inhomogeneous models 

A complete set of comparison plots from one of the inhomogeneous models tested is 
given here — a model with varying geometry/energy. This model is one in which the 
two arbitrary functions M and is take a RW form, while we vary the third function 
E. 

The correct origin parameters are i?od = 0.72, god = 0.52, which gives us a near- 
parabolic case. The extracted values are Hq « 0.71953 and qo « 0.52421. Figure [7] 
shows that the M curve plotted from our numerical output is slightly below the correct 
one; in fact there is about 2.17% error at z = 3. Although this percentage error is 
bigger than the ones we had for the homogeneous cases, this is to be expected since we 
are working with inhomogeneous data that was numerically generated. The percentage 
error is a bit larger for W , being about 26.6% as shown in Figure [H However, this 
percentage error in W is large mostly because W is quite small, and we note that the 
absolute error u\ E = {W"^ — l)/2 is about the same as before. 

From Figures [S] and [TO] we can see that the correct data and the numerical 
output are generally in good agreement for both r and is, except near the origin. 
However, this is due to insufficient accuracy in the Hq and go values deduced from 
the "observational" data at z = (3/2)i5z. Accurate values for r and ts depend on an 
accurate go value, which is particularly difficult to get at the low z values near the 
origin. A least squares estimate of origin values, using a wider range of near-origin 
data may improve accuracy here. 
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Numerics 
Correct 



Figure 7. Results of M vs. z with Hq ^ 0.71953, qo ~ 0.52421 and 5z = 0.001. 

The soUd grey curve is from the correct testing data and the dotted black curve 
is our numerical output using Runge-Kutta as the integration method. 



Numerics 
Correct 



Figure 8. Results of W vs. z with Ho ~ 0.71953, qo ~ 0.52421 and Sz = 0.001. 
The solid grey curve is from the correct testing data and the dotted black curve 
is our numerical output using Runge-Kutta as the integration method. 
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Numerics 



Figure 9. Results of t vs. z with Hq ^ 0.71953, go ~ 0.52421 and Sz = 0.001. 

The grey curve is from the correct testing data and the dotted black curve is our 
numerical output using Runge-Kutta as the integration method. 



Correct 

Numerics 

Now 



Figure 10. Results of ts vs. z with Hq ^ 0.71953, qo ~ 0.52421 and 5z = 0.001. 
The thick solid grey curve is from the correct correct testing data and the dotted 
black one is our numerical output using Runge-Kutta as the integration method. 
The solid black line is the current time (origin). 
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5. Conclusion 

We have developed a computer programme to implement the MHE algorithm. 
Given (spherically symmetric) data from standard observations for redshift, apparent 
diameter, apparent luminosity and galaxy number counts, as well as the associated 
evolution functions, true diameter, absolute luminosity and mass per source, it 
determines the metric of the (observed) universe. Its ability to reproduce the 
correct metric information has been tested via artificial data generated from both 
homogeneous and inhomogeneous models. 

We have started with a very simple case, in order to understand the key elements 
of a numerical extraction of metric information from observations. Obviously one docs 
not wish to tackle the full complexity of the problem at the start. Thus, there are still 
many improvements which can be made in both the theory and the numerical method 
used. Many considerations and effects must be included, for example, source evolution 
theories, data set completeness, different populations of sources, and more. At some 
point, a non-zero A should be considered. Also, issues like a least squares fit for the 
data near the origin in order to obtain better accuracy for i7o a-nd goj ^^nd a shorter z 
range for the series expansion in order to carry our numerics through the point Rmax , 
also need to be dealt with for the future development. Of course, a higher order 
integration method may also be needed in the future in order to sustain the accuracy 
we have so far in our numerical output out to larger z values. However, any numerical 
method we use to solve the DEs must be able to handle both known data and unknown 
functions at a discrete set of positions. Although higher order Rungc-Kutta methods 
will have a natural smoothing effect, some other form of smoothing of the data may 
be needed when we tackle real data. 

The bin size used for binning the observational data will affect the accuracy of 
the bin averages, and require attention when one works with the real data. Using 
the same bin size for the whole redshift range, will leave the higher z bins flooded 
with data, and low z bins with very sparse data. On the other hand, making nearby 
redshift bins too large may impose too much smoothing. A question for the future 
then is the optimum binning strategy, and the choice of binning versus smoothing, 
given that numerical integration is ultimately a discrete process. 

Our initial attempt at a numerical implementation of the MHE procedure has 
successfully demonstrated the viability of the basic concept, and opened the way to 
developing a more general treatment. Our current focus is on developing a workable 
numerical scheme. There are of course many relevant observational issues, such as 
luminosity functions, K corrections, different source populations, source evolution, 
bias, etc that must be incorporated in reducing the observations to the data that such 
a programme must use. These will be considered in the future. 

Current redshift surveys do not have the accuracy or completeness to enable 
meaningful metric data to be extracted at presenlQ. However, the next generation 
of surveys is expected to provide considerable improvements in accuracy and 
completeness, as well as extending to much deeper z values. Type la supernova 

tt For example, the 2dF data has large fluctuations in its number counts plot, since its thin slices were 
strongly affected by the individual clusters and voids encountered. Such large fluctuations would not 
be present in data averaged over all angles. Excessive fluctuations may cause the numerical method, 
which currently does not allow for angular variation, to break down, because we take flrst and second 
differences from the R{z) data. Small fluctuations in R will generate much exaggerated fluctuations 
in dR/dz and SR/dz'^. 
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measurements hold the promise of very good luminosity distance data in the near 
future, so the accuracy of luminosity functions and therefore number counts will be 
the limiting factor. 

Current work involves analysing the stability of the DEs, ensuring the procedure 
can handle data with statistical scatter, estimating uncertainties in the output from 
uncertainties in the observational data, and using the properties of the maximum in 
the area distance as a means to check or correct the result. 

Eventually, knowing the metric nearby will assist in analysing more distant 
observations in more than just a statistical sense, since the spacetime that the light rays 
we observe have travelled through, changes the size, brightness, frequency, position 
and shape of the images we measure. Therefore, as we probe deeper into space, the 
knowledge of the geometry of the universe around us will certainly play a crucial role 
in the data reduction of any survey in the future. With more reliable observational 
data, one may hope to achieve one of the long term objectives of the current project 
- being able to prove the homogeneity of the observable region of the universe rather 
than just assuming it in principle. 
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Appendix A. The near-Parabolic case 

One may have noticed from the three evolution equations ([l])-® that the parabolic 
evolution is actually the E ^ Q limit of the other two evolutions, which is obtained 
by writing the functions of 77 as Taylor expansions for small rj and noting that 77/%/ E 
remains finite. One can see that as we are approaching this borderline case, the 
evolution equations ([4|)-([6]) are not well-behaved numerically. Also, in reality, it is very 
difficult, if not impossible, to obtain an exactly parabolic case numerically. Hence, a 
series expansion is needed in order to have reasonable numerical results for the near- 
Parabolic case. 

Most of the series expansions for the near-Parabolic case can be found in [15]. 
However, here we will consider the hyperbolic case and give a detailed derivation 
following the approach in [15], but for obtaining the series expansion for t ~ t ~ ts 
only since this is the only one that is essential to us. Let us first introduce two new 
variables x = 2E/M^/^, and a = R/AI^''^. The parabolic limit now occurs when 
a; — > 0, while R and r remain finite. By ([4|), this requires 

77 -> and > e (A.l) 

^x 

so that the new evolution parameter e remains finite for finite r. Taylor series 
expansion expressions of r and a for the hyperbolic case using equation ([4]) are just 

(A.2) 
(A.3) 
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If we invert the scries for a by writing e in series expansion form: 

e « eo + eix + e2X^ + 632;^ H , (A. 4) 

then substituting into (jA.Sp . and solving for the coefficients e^, we get 

e - (1 - ^ + ^ a-x- - ^ a^.^ + a^x^ + ..•), (A.5) 

which we substitute into (|A.2|1 . and write it in terms of R, M and giving 

2E? ( 1 I R{2E) 3 R^{2Ef 5 R^{2Ef 



M V 3 20 M 224 AP fI52 A/3 

35 R'^{2Ef 



22528 M4 



(A.6) 



Equation (|A.6p is the r series expansion expression for the near-parabohc case. One 
can do the derivation using the eUiptic evolution equations similarly. 

Appendix B. Coefficients of the series expansions near Rmax 

Let us say that Rmax occurs at Zm, f^n{zm) = (A'^)oi a-nd we define Az = z — Zm- So 
the series expansions for p,n{z), R{z), (j){z), M{z) and W{z) have the form 

fin{z) = (/in)o + ifin)iAz + {fin)2Az'^ + {fin)3Az^ H , (B.l) 

i?(z) = R,nax + i?2Az2 + R^Az^ + R^Az^ + R^Az^ + ■■■ , (B.2) 

0(z) ^4>o + 01 Az + ^sAz^ + (t>3AzJ + • • • , (B.3) 

M(z) = Mo + MiAz + A/aAz^ + MgAz^ + • • • , (B.4) 



and 



Hence, 



and 



W{z) = Wo + WiAz + W2Az2 + • • • . (B.5) 



dR 

—— = 2i?2Az + 3i?,3Az^ + 4i?4Az3 + Si^gAz"* + • • • , (B.6) 
az 

d^R 



2R2 + 6R3AZ + URiAz^ + 20i?5Az3 + . . . . (B.7) 
az^ 

And the expressions for the coefficients of the series expansion for Af (z), 0(z) and 
W{z) are given below: 

Mo = , (B.8) 

Ml = Ml , (B.9) 
f^_^S^\M^_R^_8fi^ 

V 1 + ^™ (A«)o J 2 2 ^ ' 
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(B.ll) 
(B.12) 
(B.13) 
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Appendix C. The near origin expressions 



(B.18) 



We do series expansion of equations (A.l) and (A. 2) in Appendix A of MHE near the 
origin, obtaining 



1 (3 + go) 2 4 + go + go2 ^ 
z — z H" z 



Ho 



2Ho 



2Ho 



and 

47r/in ~ z'^ - + ^3 ^ 3go(15go^ + 14go + 13) ^4 ^ 



4i?o 
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We test the accuracy of the generated Hq and go values using z ~ {?)/2)5z^ and 
R and Anjln values at this same z since this is where we have the first complete set 
of data according to available observational data and the way we discretisc our DEs, 
and therefore, find the combination with the most consistent accuracy for different go 
values. We find that in general, using both equations (|C.1|) and (|C.2[) with the same 
number of terms gives us better accuracy. We can then get expressions for i?o and go 
near the origin in terms of z, i? and An fin only. The results are 



and 



J (47r/i7i)2 + 36i?2(2z - 1)2 + 12i?(47rA7i)(10z - 7) 

go = ^ 

24zi? 

+ 4^An + 6i^(l-2z) 
24zi? 

3goz^(l-2goz-2z) 
= ■ (^-^^ 

Now we have a way of determining the origin values for Hq and go from the 
data. If we need to generate values for r, 0, M and W dX z = {l/2)5z and the origin 
numerically from the RW expressions given in Appendix A in MHE, then one can 
perform series expansions of them too, since the values of z are small. They are: 

90 3 3go(H-go) 4 3go(3 + 2go + Sgg^) 5 

Ivl fiy/ ~ z — z H z 

Ho 2 Ho AHq 

go(28go3 + I2q^^ + 15 go + 25) ^ . . . . (^.5) 



8i/o 

I R 3 , R„ 2 I R„_ I fi ^ 



_ 2 + go W + 4go + 6 5go=^ + 6go^ + 6go + 8 _^ 
ORW FT 77 ^ ;TI7 ^ ^ : (.'-^■Dj 



Ho Hq 2Ho 2Ho 

2Erw « (1 - 2go)z2 - go(l - 2qo)z^ + " 2go)(5go2 - 2go + 1)^^ 

-i(l-2go)(7go3-4go2 + l)z5 + ... ; (C.7) 
1 2 + go 2 , 3go' + 4go + 6 3 Sgo''^ + 6go' + 6go + 8 . 

^«^"75o'-^i^^ + — Wo — ' Wo ' +---;(c-8) 

and these expressions are valid for any go. For go < 1/2 and go > 1/2 respectively, 
the T series expansions for the RW ec^uations take the form: 



^/^^ + g°ln^T7T=i^J 1 2 + go , 

ffo(l-2go)3/2 Ho 2Ho " 
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6ffo 8iIo 
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